# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <cmath>

using namespace std;

#include "lagrange_1d.hh"
#include "fe_evaluators.hh"
#include "utils.hh"

int main ( )
{
  typedef double BaseType ;
  typedef fadbad::F<double> GradType ;
  typedef fadbad::B<fadbad::F<double> > HessType ;

  typedef GradientEvaluator<QuarticLagrange1D, BaseType> GradEvaluator;
  typedef HessianEvaluator<QuarticLagrange1D, BaseType> HessEvaluator;

  timestamp ( );

  cout << "\n  Test the FEM evaluation functions.\n" ;
  
  int test_n = 10;
  Point<SPATIAL_DIMENSION, BaseType> pt ;

  cout << "\n  The Basis should sum to 1 and X derivatives should sum to 0.\n\n";
  cout << "  \t\t(Phi_i, dPhi_i/dX_j) ... \t\t(sum(Phi), sum(dPhi/dX)):\n\n";

  GradEvaluator gradeval;
  
  Point<SPATIAL_DIMENSION, BaseType> grad ;
  for (unsigned int j = 0; j < test_n; j++ )
  {
    pt[0] = generate_random();

    BaseType basis_val = 0.0, sum_x = 0.0, sum_dx = 0.0;
    cout << "  ";
    for (unsigned int i = 0; i < FEMTraits<HessEvaluator::FiniteElementType>::N_Nodes; i++ )
    {
      gradeval(i, pt, basis_val, grad);
      cout << setw(8) << "(" << basis_val << ", " << setw(8) << grad[0] << ")" ;
      sum_x += basis_val;
      sum_dx += grad[0];
    }
    cout << setw(8) << "(" << sum_x << ", " << setw(8) << sum_dx << ")" << endl ;
  }

  cout << "\n  The Basis should sum to 1 and X derivatives and double-derivatives should sum to 0.\n\n";
  cout << "  \t\t(Phi_i, dPhi_i/dX_j, d2Phi/dX_ij) ... \t\t(sum(Phi), sum(dPhi/dX), sum(d2Phi/dX_ij)):\n\n";

  HessEvaluator hesseval;
  
  PointList<SPATIAL_DIMENSION, BaseType> hess ;
  for (unsigned int j = 0; j < test_n; j++ )
  {
    pt[0] = generate_random();

    BaseType basis_val = 0.0, sum_x = 0.0, sum_dx = 0.0, sum_ddx = 0.0;
    cout << "  ";
    for (unsigned int i = 0; i < FEMTraits<HessEvaluator::FiniteElementType>::N_Nodes; i++ )
    {
      hesseval(i, pt, basis_val, grad, hess);
      cout << setw(8) << "(" << basis_val << ", " << setw(8) << grad[0] << ", " << setw(8) << hess.loc[0] << ")" ;
      sum_x += basis_val;
      sum_dx += grad[0];
      sum_ddx += hess.loc[0];
    }
    cout << setw(8) << "(" << sum_x << ", " << setw(8) << sum_dx << ", " << setw(8) << sum_ddx << ")" << endl ;
  }


//
//  Terminate.
//
  cout << "\n  Normal end of execution.\n\n";

  timestamp ( );

  return 0;
}
